Exact duality and dual Monte-Carlo simulation for the Bosonic Hubbard model 
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We derive the exact dual to the Bosonic Hubbard model. The dual variables take the form of 
conserved current loops (local and global) . Previously this has been done only for the very soft core 
model at very high density. No such approximations are made here. In particular, the dual of the 
hard core model is shown to have a very simple form which is then used to construct an efficient 
Monte Carlo algorithm which is quite similar to the World Line algorithm but with some important 
differences. For example, with this algorithm we can measure easily the correlation function of 
the order parameter (Green function), a quantity which is extremely difficult to measure with the 
standard World Line algorithm. We demonstrate the algorithm for the one and two dimensional 
hardcore Bosonic Hubbard models. We present new results especially for the Green function and 
zero mode filling fraction in the two dimensional hardcore model. 

PACS numbers: 75.fONr, 05.30 Jp, 67.40.Yv, 74.60. Ge 



I. INTRODUCTION 

Since its discovery by Kramers and WannierS, dual- 
ity has continued to play an important role in physics. 
Originally, duality related the "high temperature" phase 
of one system to the "low temperature" phase of another. 
But it has long been realized that duality is more general 
than what is encountered in thermal statistical physics. 
In the case of quantum phase transitions, for example, 
high and low temperatures are replaced by high and low 
coupling constants which determine the size of quantum 
fluctuations. As an example in the thermal case, the two 
dimensional Ising model is self-dual i.e. its low temper- 
ature phase is mapped onto the high temperature phase 
of its dual which is also an Ising model. This allowed 
Kramers and Wanniero to determine its exact critical 
temperature. The nature of the dual and the question 
of self duality depend on several factors like the symme- 
try group of the action, the Jtype of model (spin, gauge 
etc) and thcj-||dimensionality.D Duality has found many 
applications.crB For example, a very interesting recent ap- 
plication was to show that the bosonization of fermionic 
systems in any number of dimeasions can be effected us- 
ing the duality transformation.Q 

In quantum systems, such as the Bosonic Hubbard 
model, an approximate duality transformation has been 
used to express the partitioa function in terms of con- 
served electric current loops.B This transformation is ap- 
proximate in two ways: (1) The boson density is assumed 
to be very high, in other words the model is necessaxily 
very soft core, and (2) the Villain form of the actionH is 
assumed. Because of this, the dual model thus obtained 
is not even the dual of the Hubbard model but that of 
the Quantum Phase Model (QPM) . Still, this formulation 
was used very productively to study various aspects of the. 
ground state phftS« diagram of the QPM numericallyEI'Ll 
and analytically.^ However, to study the bosonic Hub- 
bard model with no approximations (other than Trot- 
ter discretization), one still had to use one of the other 
available algorithms. The World Line algorithnJlJ proved 



very usefuj. in one and two dimensions with or without 
disorder .IlirEj Improvements on this (in the hard c0|Cej 
case) came in the form of the loop (or cluster|)_algorithinll^ 
and the continuous time cluster algorithmO These two 
algorithms were designed to eliminate long relaxation 
times (critical slowing down) at half filling and in the ab- 
sence of longer range interactions (like next near neigh- 
bour). When these conditions are not met these algo- 
rithms become extremely inefficient. In addition to ques- 
tions of inefficiency, a disadvantage shared by these al- 
gorithms is that even though one can calculate, rather 
easily, the superfluid density and other phvsiyal Quanti- 
ties such as the energy and compressibilityDu'EJIlj it is 
very difficult, if not impossible, to measure the correla- 
tion function of the order parameter (the Green function, 
G) whichpis a very interesting quantity. The "worm" 
algorithnJl3 can be used to calculate G but seems to 
converge very slowly in two dimensions. The--|recently 
developed stochastic series expansion methodic is very 
promising. It overcomes the inefficiencies of the loop al- 
gorithm, but whether it is efficient for measuring G is to 
be seen. 

It is sometimes desirable to do the simulation in the 
canonical ensemble especially when studying questions of 
phase separation. To this end, we will first show here how 
to implement the dualitv^transformation exactly following 
the method of reference H. We will then construct a QMC 
algorithm valid in the canonical and grand canonical en- 
sembles and with or without disorder. The efficiency of 
this algorithm does not deteriorate when the system is 
doped away from half filling or in the presence of longer 
range interactions. It also has the advantage of being 
able to calculate, rather easily, G in the grand canonical 
and canonical ensembles. 

The paper is organized as follows. In section II we 
show the duality tranformation leaving the details for the 
appendix. In section III we simplify the results and con- 
struct the QMC algorithm for hardcore Bosons. In sec- 
tion IV we discuss some algorithmic improvements and 
details and show tests of the algorithm in the one dimen- 
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sional bosonic Hubbard model in the hardcore hmit. In 
section V we present new simulation results for the two 
dimensional hardcore Bosonic Hubbard model with near 
and next near neighbors. Section VI is for conclusions. 



II. DUALITY TRANSFORMATION OF THE 
BOSONIC HUBBARD MODEL 

A. The Bosonic Hubbard model 

We are interested in simulating models whose Hamil- 
tonian has the form 

H = -tJ2 + alar) + VoJ2 ^ ^ (1) 

(rr') r 

+ Vi flrhr' + V2 flrflr' + • • • — /i^^Tlr. 

(rr') {{rr')) r 

In this equation, t is the transfer integral (the hopping 
parameter) which sets the energy scale (in the simula- 
tions wc set it equal to 1), fi is the chemical potential, 
Vo, Vi, and V2 are respectively the contact, the near and 
the next near neighbor interactions. The dots stands 
for other longer range interactions that can be added. 
On the square lattice, the next near neighbor is chosen 
along the diagonal, while in one dimension the choice is 
the obvious one. a^- and al are destruction and creation 
operators on site r satisfying the usual softcore boson 
commutation laws, [ar,aj/] = Sry , and fir = alar is the 
number operator on site r. (rr') and {{rr')) label sums 
over near and next near neighbors. The first term which 
describes the hopping of the bosons among near neighbor 
sites gives the kinetic energy. 

The quantum statistical mechanics of this system is 
given by the partition function Z, 



0H], 



(2) 



where /3 is the inverse temperature and the sum runs 
over a complete set of states. The expectation value of a 
quantum operator is given by 



{O) = ^TiiOe-f'"). 



(3) 



B. Path integral representation 

We now need to find a c-number representation of the 
partition function, which can then be incorporated into 
a simulation algorithm. There are several such represen- 
tations. 



The occupation number representation is defined by 
the eigenstates |n) of the number operator alar, 



|n) = 



n 



(4) 



|0), 



(4) 



where |n) stands for a complete set of n^, |n) = 
\ni,n2, • ■ ■) where rir is the number of bosons at site r. 
In this representation the resolution of unity reads 



(5) 



Another common representation uses coherent states 
1$), i.e. eigenstates of the destruction operator. 



($14 = ($10*, 



(6) 



where the eigenvalues are complex numbers defined on 
the sites of the lattice. We completely define this state 
with the vector $ = (cii, . . . , (f>r, . . .). 

In terms of the occupation number representation vac- 
uum |0) the coherent state |$) is defined as 

|$)=exp|^^</..atj|0), (7) 

so that its projection on an occupation number state is 



With this normalization we obtain the inner product of 
two coherent states 



(*!$)= exp||^VrVrj , 



(9) 



and the resolution of unity 

1= /n^«-p(-Ei'^^i')i^)<^i- (10) 

r \ r / 

The interaction and chemical potential part of the 
hamiltonian are more easily expressed in the occupation 
number representation while a coherent state representa- 
tion is well suited to the treatment of the kinetic term. 
We will use both successively. We write the partition 
function, Eq. 2, as 



{n} 

{n} 



-ISHi 



n) 



-SH 



n), 



(11) 



where 6 = 13/ L^, and is a large enough integer such 
that 5^1. We express the partition function this way 
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because we can now express the exponentials in Eq. 11 
in a form suitable for easy evaluation. Between each pair 
of exponentials introduce the resolution of unity, Eq. 5. 
We find 



{n,} 



As 5 is small, 



^~&H ^ g-5(T+V) 



(n3|e-*^|n2)(n2|e-^^|ni) 



(12) 



g-5y/2g-5Tg-5y/2^(^(^3)^ (13) 



where T is the kinetic energy and V contains the inter- 
action and chemical potential terms. 
The partition function now becomes 



-SV{iit) 



(ni|e-*^|nz.J 



(14) 



The error introduced by neglecting commutators of T 
and V is of order Lt5^ = /3t^ for the partition function. 

We now introduce unity operators expressed in terms 
of coherent states (Eq. 10) around the kinetic terms. 

^ = E /n f ^!^^4!^e-(i'^'-.-i' + i^'-.-i')) 

X n (ni||.i)($i|e-*^|$i^)($ijni,) 

T 

X (n2|l'2)(l'2|e-'^|$i)($i|ni). (15) 

Notice that we need to introduce two sets of coher- 
ent states, the $ and $. We can avoid introducing two 
different representations (occupation and coherent) and 
only use coherent states. But then, we would have to de- 
couple the interaction terms in order to integrate them 
out. At the end this would lead to the same results, but 
it would be more difficult to recover a simple expression 
for Z. 

At this stage, it is customary to neglect terms of order 
2 and higher in 5 in the expectation value of exp(— (5T) 
since their contribution goes to zero with 6. While this 
is correct and would still allow us to calculate thermody- 
namic quantities such as {E)^ it would not be sufficient 
to calculate the Green function, G{\r — r'\) = {ara\i), at 
distances |r — r'| > 1 (see below). 

For this reason, it is very important to keep higher 
orders. We start by keeping all such higher order terms 
of exp(— (5r). Since T is quadratic 



T ■ 



jtKE 



(16) 



where K-js a c-number matrix, we can use the following 
identityO 



= {^\ : exp I 



-at fl - e-*K] a) 



(17) 



where : O : indicates that the operatoij-C* is normal or- 
dered. Eq. 17 may now be expressed asliZI 



(*|e- 



(Sa^Kal 



= (^'1$) exp(-*t [1 
= e*'* exp [1 _ e-^K] (^g) 

Expanding the exponential term in Eq. 18, one finds 
5^K^ 



e exp 



SK 



2! 



3! 



4! 



(19) 



For simplicity, we will illustrate how to proceed by us- 
ing one dimensional L^j-site lattice with periodic bound- 
ary conditions. All nearest neighbors {rr') pairs are la- 
beled by choosing r' = r + 1 and summing over all r. The 
same can be easily done in any dimensionality. We write 
K as K = — tKi where Ki is the matrix that connects 
nearest neighbors: If r and r' are nearest neighbors then 
Kl^ = 1, otherwise K^^ =0. It is easy to show that 

K2 = (-i)2(2I + K2), 

= (-i)4(6I + 4K2 + K4)--- (20) 

where generates hops of L lattice spacings. 

Introducing this into Eq. 19 and keeping only the lead- 
ing order term in S for each matrix, K^, we get 



exp 



vl/t 



2! 
3! 



(21) 



4! 



Substituting the preceding expression into Eq. 15, us- 
ing Eq. 8 to express the scalar products, and neglecting 
higher order terms, we get the following expression for 
the partition function 



Z = 



n 



r,r 



7'($,$,n), 



(22) 



where the weight V is given by 



7'(a>,|.,n)=[][e 



(23) 




2! 



3 



We reiterate that as i5t ^ 0, the contributions of mul- 
tiple hops, i.e. K.L where L > 1, to the energy and 
other local physical quantities vanish since they are of 
higher order in St. However, this is not true for non lo- 
cal quantities such as G(r) — {aoal) for r > 1. Keeping 
the higher order terms will give us the generating func- 
tional for this correlation function. Notice, however, that 
we have kept only the terms necessary to calculate each 
quantity to leading order in S. Therefore the algorithm 
will have systematic Trotter errors of order (3S. 



C. Dual formulation 

In order to perform the integration over the original 
variables and obtain the dual formulation of Z, we express 
the (j)r,T and (pr^r as, 

0r,r = /Or,r exp(i6'r,r), (24) 

which gives the following expression for Z 
1 



Z 



',T ^r,T ) "'r,T 



X exp [ Pr,r + lPr.re~'^^'-'^+^ ~ ^''^^^ 

X exp f 5tpr,^+i/9r+i^^e~*(^''^^+i ~ ^r+\,T 



(25) 



X exp 



X exp 



2! 



-Pr+2,T + lPr,Te 



2! 



Dots stand for all longer range hopping terms and J T)"^ 
stands for 



E / n 



\PrA^ + |Pr,rP) 



-Er^^(nr). (26) 



As expected, only gauge invariant phase differences ap- 
pear in Eq. 25. Therefore, instead of performing the inte- 
grals over the gauge dependent site phases, 9r,T and 9r,T^ 
we express the partition function in term of the following 
gauge invariant quantities 



e° = 

= 



'r.T ^r,Ti 



(27) 



We use indeces and 1 to label the imaginary-time and 
X directions. In higher dimension, we would have 2 cor- 
responding to the y direction, 3 to z... 

This variable change results in a non-triviaL Jacobian 
which was derived in reference 3. It was shownci that this 
Jacobian is only the product of two kinds of constraints 
on the gauge variables and that it can be written as 



^= E E n 



(28) 



where the integer valued variables, A^^^ and J/:*^, rep- 
resent conserved currents and are, in fact, the dual 
variableg3. The origin of these constraints is that the 
sum of the bond variables along any directed closed path 
must be zero as is seen from their definition in Eq. 27. 
0^ does not appear in the hamiltonian and so has been 
eliminated in the Jacobian. For a detailed discussion of 
these constraints and their relation to duality, the strong 
coupling expansion etc see Ref 3. 

The global currents, A^^, traverse the system from 
one end to the other and are required by the periodic 
boundary conditions. A'^ is a global current in the time 
direction and, therefore, does not depend on the coordi- 
nate r. Similarly, A^^, the global current in the x direc- 
tion depends only on r (in two dimensions it would also 
depend on y). The local currents, Jl^^., form topologi- 
cally trivial closed loops, i.e. they do not wrap around 
the system from one end to the other. The sum of the 
J^^ configurations is restricted to be only over conserved 
current configurations: The total current entering a site 
equals the current leaving. This is indicated by the prime 
on the sum. We will see below that this current has a 
simple physical interpretation. 

With this variable change and Jacobian, Z becomes 



Z = 



\Pr,T Pr.T j ' ^ ' ' 



-zeo,) 

X exp ( Stpr+l,T+lPr,re 



X exp \ Pr,r+lPr,Te ~ " r,T , ^ (29) 



X exp 5tpr,T+lPr+l,TG 



X exp 



X exp 



2! 



-Pr+2,T+lPr,: 



"*(Qr-(-2,r + ©r-|-l,T + ©r.r) 



2! 



-Pr,T+lPr+2,T 
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where 



/ 



(30) 



n 



X e 



Pr.T Pr.T 



-(IPr.rP + IPnrH 



Integrating the original variables, 6, p, and n leaves 

only the integer valued dnal variables. The details are 
shown in the appendix. The result of this integration is 
a complicated expression for the dual partition function 
involving several integer valued fields and constraints. 

One can greatly simplify this expression by explicitly 
solving the constraints appearing in it (these constraints 
come from the integration of 0), that is, by finding the 
allowed configurations of Nl^!-^^ and J^^ and their asso- 
ciated Boltzmann weight. In the appendix we illustrate 
this for softcore and hardcore bosons but concentrate on 
the hard core case where the dual partition function takes 
a particularly simple form. 



III. QUANTUM MONTE CARLO ALGORITHM 
IN THE HARD CORE LIMIT 

A. Partition function in the hardcore hmit 

As shown in the appendix, we can express the par- 
tition function and other physical quantities in terms 
of only one field, the total conserved hard-core current 
N!^^ + J^^. By hard-core, we mean that the total current, 
TV^.^ + J^^, traversing a bond in the imaginary time di- 
rection takes the values or 1, while in a spatial direction 
it takes the values 0, ±1. We have also chosen, without 
any loss of generality, to keep ^ J^^ = 0. Notice that 
the total time currents are never negative: The bosons 
cannot go backward in the imaginary time direction. 

The partition function in the grand canonical case 
reads 



Z= ^ ViN,J), 

{J,N} 

Y(<5^)^^^^ 




(31) 



where N]^ is the number of jumps of length L on the 
lattice, the number of near neighbor time currents, 
and the number of next near neighbor time currents. 



Time 



Space 

FIG. 1. Example of a configuration of the conserved Boson 
current on a Id lattice. 



The chemical potential appears in an exponential of 
l3J2rT ^rr- number of particles in the system is 

simply 



(32) 



which demonstrates that the current simply represents 
the current of bosons crossing the lattice in the positive 
imaginary time direction. So we have current lines cross- 
ing the lattice and the sampling of the configurations is 
achieved by deforming these lines. If one wants to work in 
the canonical ensemble, one simply has to fix the number 
of current lines. 

The effect of interactions is simply given by multiply- 
ing the weight by e~^^^ (e~''^2 ) for every pair of near 
neighbor (next near neighbor) time currents. 

The term containing t, coming from the expansion of 
the hopping term, also has a simple physical interpre- 
tation: Kinetic energy arises when there are jumps, i.e. 
spatial currents. This term can also be written as 



n 



Nr 



(33) 



The exponent Nh of 5t is the number of non zero spatial 
currents Nh = = J2r t l-^r + Jr,r\ t>ut there is a 

combinatorial prefactor arising from multiple jumps. 

We therefore have a simple physically intuitive expres- 
sion for Z which allows easy evaluation of the Boltzmann 
weight. For example in Fig. 1, there are 7 single jumps, 
one double jump, and one triple jump. There are also 
one pair of near neighbor time slices and one pair of next 
near neighbor. So the weight associated with this config- 
uration is 



(St) 



12 



2! 3! 



(34) 
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B. Physical quantities 

As shown above, we have for the particle number, n^.r, 

{nr,r) = + Jr,r) ■ (35) 

Then the density-density correlation function reads 



ill 



(n(ri)n(r2)) = 



ri,T / 



r2,T/ 



(36)L 



and the structure factor, S{k), is its Fourier transform. 

We now focus on the canonical case. The energy is 
given by (E) — — ^ In Z and therefore by 



(E) - - (SViNv, + 5V2NV, - Nh) 



(37) 



The meaning of the three terms is transparent. The first 
two are the near and next near neighbor potential ener- 
gies, and the last, which comes from the jumps, is the 
kinetic energy. 

The superfluid density is related to the winding num- 
ber, W, byBlia 



Ps 



2t(3 ■ 



(38) 



where W, given by J^t -^t j changes with the global spa- 
tial currents N^. 

Finally, we are able to calculate easily the Green func- 
tion, i.e. the correlation function of the order parameter. 



G{L) 



1 



which, in isotropic cases, equals 



(ar^ral^^^ + ar+L.ralr) ) ■ (39) 



G{L) = 



1 



LtLx 



(40) 



This can be calculated directly from the generating func- 
tional, Eq. 31, 



G{L) 



1 



d 



'2i L-fLj^ 



d 



5H^ 



■In Z, 



(41) 



which can be easily shown to equal 
1 LI 



G{L) = 



-Nr 



(42) 



We now see the importance of keeping multiple jump 
terms: The only configurations that contribute to G{L) 
are those that contain jumps of length L. This is physi- 
cally satisfying: There are correlations between two sites 
if they are directly connected by a jump. This quan- 
tity is also simple to calculate, we just need to count the 
number of multiple jumps on the lattice. 



(a) (b) 

FIG. 2. The two kinds of trial moves for the one dimen- 
sional model: Local moves (a) and global moves (b). 

IV. ALGORITHMIC DETAILS 

A. Evolution of the system 

We begin by choosing an initial configuration of 
Bosons. This is done by introducing a current line for 
each particle i.e. by putting the corresponding global 
time currents to 1. Then the simulation proceeds by 
modifying these current lines using a Metropolis sam- 
pling scheme. 

In one dimension, we can explore the whole phase space 
by simply applying two kinds of trial moves. First, we 
introduce local moves (Fig. 2(a)): We visit sites sequen- 
tially and try to add a small (clockwise or anticlock- 
wise) current loop. If this does not introduce currents 
greater than 1 or negative time currents, we accept or 
reject probalistically according to the Boltzmann weight, 
Eq. 31. In addition to these local moves, we have global 
ones (Fig. 2(b)) where we change the value of a spa- 
tial global current traversing the system. These moves 
change the winding number. The global moves have very 
low acceptance probability, which means that in practice, 
one cannot efficiently sample all the winding number sec- 
tors. 

In order to do simulations in the grand canonical en- 
semble, one can introduce global moves along the time 
direction which will add or remove particles. One can 
also use the chemical potential term to introduce a dis- 
ordered external potential by imposing a site dependent 
fi. 

In two dimensions, we need to introduce two more 
types of trial moves to get topologically different con- 
figurations and ensure ergodicity. These moves are the 
one particle twist and two particle exchange shown in 
Fig. 3. 
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B. Reweighting 





(a) (b) 
FIG. 3. Additional trial moves for the two dimensional 
case: One-particle twist (a) and two-particle exchange (b). 



(a) (b) 
FIG. 4. A configuration of the conserved current for one 
boson (a) with four single jumps and (b) with one double 
jump and two single ones. 



As we have seen, multiple hopping terms are necessary 

for the calculation of the Green function, G. Consider 
the two configurations of Fig. 4, (a) contains only single 
jumps and has a weight {St^, (b) contains one double 
jump and two single ones and so its weight is {5tY/{2\). 
Both of them give a contribution of the same order in 
5 to Z. In the first case, since the position of each of 
the four jumps is arbitrary, the total number of such 
configurations grows as V^. Consequently, they give a 
finite contribution of order L'^{St)'^ = (pt)^. For (b), 
only three jumps have arbitrary positions, which gives a 
final contribution of (5^)4 = lpt)^St. Taking the ^ ^ 
limit, we see that the configurations with multiple jumps 
give no contribution to the partition function. 

For the same reason, all configurations with jumps of 
size L + I and many jumps of size L, do not give a finite 
contribution to G{L). The only important configurations 
for G{L) are those with one jump of size L and L single 
jumps. 

Sampling the configuration with the Boltzmann weight 
(Eq. 31), the probability of having multiple jumps falls 
to zero as S gets smaller. This causes a problem: We 
want small 6 in order to reduce the Trotter discretization 
error which is proportional to S. But in that limit the 
number of multiple jumps falls exponentially making it 
extremely difficult to calculate the Green function G. 

To overcome this problem and calculate Green func- 
tion efficiently, we introduce a reweighting scheme. We 
replace the weight P by a weight V-ji designed to enhance 
configurations with multiple jumps. 

Equations 31 and 42 show that a natural choice is to 
keep Vii = V if there is no jump. When there are jumps 
we take 



Vr 



[{6t)M-i/M !] 



(43) 



where M. is the length of the biggest jump on the lattice. 
Since this changes the Boltzmann weight, we should com- 
pensate when measuring physical quantities. For the en- 
ergy, we get (the sums run over all allowed configurations 
of the current) 



{E) = 



1Y.{5ViNv,-Nh)V 
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(44) 



1 E mNv, - Nh) ^Vr 



1 E {SViNy, - Nh) [{5t)^-^/M !] Vr 

(3 Em^-'/M\]pR 

1 {{SViNv, - Nh) [iSt)M-^/M !]) 



Vr 



m)M-yM\])T,^ 

and for the correlation function 
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FIG. 5. Number of jumps of different lengths with and 
without reweighting for a system of 16 bosons on 32 sites 
at /3 = 4. 



G{L) 



L\ {[{5t)^-^IM^]NL) 



Vr 



(45) 



Vr 



This new reweightcd algorithm increases dramatically 
the efSciency of the calculation of G. Figure 5 shows the 
number of jumps as a function of the length of the jump 
using the original and rcweighted algorithms. As one can 
see, it falls exponentially with length in both cases. But, 
for the reweighted algorithm the slope is much weaker 
allowing frequent long jumps. 

Reweighting favors long jumps which we can exploit to 
change the global spatial currents and thus the winding 
number. This algorithm then gives us exact results (after 
extrapolation to 5 = 0) even for small size systems where 
the effect of non-zero W cannot be neglected. 

That does not mean that using the Boltzmann weight 
without reweighting is of no interest. For big systems, 
non-zero winding numbers will have very small effect on 
local quantities and if one only wants these, there is no 
need for reweighting. 



C. Tests of the algorithm 

In order to test the algorithm, we compare its re- 
sults with those of exact diagonalization on a small one- 
dimensional system (4 bosons on an 8 sites lattice at 
/3 = 4 with only hard core interactions). We extrapolated 
our results to the 5 = Q limit before comparing (Fig. 6). 
The systematic errors are, as expected, clearly linear in 5. 
The extrapolated results are in excellent agreement (less 
than one percent error) with exact ones. Fig. 7 shows 
comparison for density-density correlation function and 
G. 
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FIG. 6. Extrapolation to 5 = of the Green function at 
different distances for a system of 4 bosons on 8 sites at /3 = 4 
with only hardcore interactions. 
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FIG. 7. Comparison of Green function and density-density 
correlation functions obtained with the dual algorithm and ex- 
act diagonalisation techniques. The system contains 4 bosons 
on 8 sites with (3 = 4 and only hardcore interactions. 
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V. APPLICATIONS 



As discussed above, with this algorithm (as with oth- 
ers before it) we can measure various physical quanti- 
ties, such as the energy, the superfluid density, ps, the 
density-density correlation function and its Fourier trans- 
form, the structure factor. What is new here is the 
ease with which we can measure the Green function, 
{G{L) = a^{L)a{0)) which, among other things, gives 
the number of condensed particles (occupation of the zero 
momentum mode) . If one defines 



G(0) 



(46) 



then the Fourier transform of G{L) equals the density of 
particles in the kth mode. In one dimension 



1 



^e'''^G{L) = N{k) = 



(alak) 



where 



(47) 



(48) 




FIG. 8. Green function and density-density correlation 

function (inset) versus distance. 72 bosons on a 12 x 12 lattice 
with Vi = 0, V2 = 3 and /? = 6. The bosons form a striped 
incompressible solid. 



Using G{L) and A^(0), one can easily obtain valuable in- 
formation about the presence (or absence) of off-diagonal 
long range order (ODLRO), i.e., phase coherence at long 
distances. 

We studied with this algorithm a two-dimensional 

12x12 system of hardcore bosons at various fillings and 
interactions. Wc took /3 = 6, which is a low enough 
temperatm'c to access ground state properties. How- 
ever, although the information is present in the configu- 
rations generated by this algorithm, we did not measure 
the Green function off-axis. We only measured it along 
the two lattice axes to reduce measurement time. 

In Fig. 8 we show results for a 12X12 system at half fill- 
ing and strong next near neighbor repulsion, = 0, V2 = 
3. The nun repulsion makes the presence of next near 
neighbours too costly. Consequently, the Bosons orga- 
nize themselves into stripes with a (7r,0) (or (0, tt)) or- 
dering vector for the structure factor. In other words, the 
Bosons form an incompressible gapped solid with stripes 
along the x ov y axis. This is shown very clearly in the 
density-density correlation function shown in the inset of 
Fig. 8. Also shown in Fig. 8 is the Green function along 
the X and y axes. Notice that G{L) = {a{Q)a^{L)) falls 
faster along the stripes (x-axis) than in the transverse 
direction (y-axis). This is because the energy cost for 
transverse hops, while high, is still finite. On the other 
hand, the cost for hops along the stripes is infinite be- 
cause such hops would always lead to multiple occupancy 
of sites. In both cases, G{L) goes to zero at long distance 
indicating the absence of superfluidity. This is confirmed 
with direct measurement of Ps = {W'^)/2tp. 
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FIG. 9. Green function and density-density correlation 
function (inset) in the superfluid phase (72 bosons on a 12 x 12 
lattice with Vi = 0, V2 = 1 and /3 = 6). The system is 
isotropic and superfluid. 
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FIG. 10. Green function and density-density correlation 
function (inset) in the supersolid phase (68 bosons on a 12 x 12 
lattice with Vi — 0,V2 = 3 and l3 — 6). The system shows at 
the same time superfluidity and long range density modula- 
tions. 



Reducing the nnn repulsion allows the long range or- 
der to melt and be replaced by ODLRO. This is shown 
in Fig. 9. The inset shows the density-density correla- 
tion functions along the x and y directions to be equal. 
The system is therefore spatially isotropic. This is con- 
firmed by the isotropy (withing error bars) of G{L) along 
X and y. Also, notice that G( ^Xis finite at long distance 
indicating symmetry breakinglla (nonzero (a)) and thus 
nonzero superfluid density which we also confirmed by di- 
rect measurement of ps — 0.304± 0.017 from the winding 
number fluctuations. 

If, instead of reducing V2 at half filling, we now repeat 
the simulations that gave Fig. 8 but doped away from 
half filling, p = 0.472, we observe more structure devel- 
oping in the Green function. This is shown in Fig. 10. 
As before, the inset shows the density-density correlation 
function along x and y. As in Fig. 8, this shows clearly 
the presence of stripe order along the x-axis. However 
we also see that G{L) does not decay to zero at long 
distance. Instead, along the x-axis, G{L) resembles that 
in Fig. 9, i.e. tending towards a nonzero value at long 
distance indicating ODLRO and superfiuidity. This is 
again confirmed by measuring directly the fluctuations 
of the winding number giving p^ = 0.027 ± 0.003. Along 
the y-axis, G{L) shows an oscillatory structure: At long 
distances, G{L — even) is greater than G{L = odd) and 
is almost equal to G{L) along the x-axis. This shows 
clearly that there is superfiuidity along the y-axis, trans- 
verse to the direction of the stripes. Direct measure- 
ment gives p| — 0.015 ± 0.005 which is less than as 
may have been expected. The fact that both, long range 
density wave structure (diagonal order) and superfluidity 
(ODLRO) exist at the same time shows this phase to be 
a striped supersolid. Thiii phase has been shown to hasffi 
a finite critical velocitytj and not to phase separateEj. 
It is therefore stable. Doping above half filling will yield 
similar results: This model has particlc-hole symmetry 



around p = 1/2. 

The wave structure in G{L) along the y-axis is easy 
to interpret physically. When the system was at half 
filling, the Bosons could not hop along the stripes (since 
no double occupancy is allowed) and could only hop with 
difficulty transverse to the stripes into the empty sites 
due to the high energy cost. When the system is doped 
below half filling, the Bosons can easily hop along the 
stripes into neighboring holes. They still cannot easily 
hop an odd number of lattice spacings into the mostly 
empty stripes due to the high cost of nnn interactions. 
However, by hopping an even number of lattice spacings, 
they can fill in holes in the Boson stripes, thus aquiring 
near neighbors at zero potential energy cost and avoiding 
next near neighbors which cost a lot of potential energy. 
This explains the smaller values of G{L) for odd L along 
the y-axis. 

This also brings out an important point. The inter- 
pretation of a striped supersolid is not simply that the 
stripes form channels along which the Bosons delocalize 
and the superfiuid flows. While this does indeed happen, 
Fig. 10 shows clearly that there is also delocalization and 
superfiuidity in the direction transverse to the stripes. 
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FIG. 11a. The filling fraction of the kx = Q mode versus 
the interaction Vi (V2 = 0) for different sizes. The systems 
are half filled and /3 = 6. 



.er the case 7^ 0, V2 = 0. As is well 
for Vi > 2, p = 0.5 the Bosons are or- 



We now 
establishedEini 
ganized into an incompressible checkerboard solid, while 
for Vi < 2 they delocalize into a superfluid. In other 
words, for large Vi, the Green function decays rapidly 
to zero as in Fig. 8 (but is isotropic in x and y), while 
for small Vi it is nonzero for large separations signalling 
spontaneous symmetry breaking. 

Measuring G{L) as a function of Vi enables us to calcu- 
late the filling fraction of the modes, {a'^{kx, ky)a{kx, ky)) 
(Eq. 47), as a function of interaction. Strictly speaking, 
this requires G{\x — x'\, \y — y'|), i.e. off-axis as well as 
on-axis contributions to G. But, as explained above, we 
measured only on-axis contributions to G. We use this to 
obtain an estimate of the filling fraction as follows: Since 
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FIG. lib. The filling fractions of Fig. 11a extrapolated to 
infinite systems. Inset: Filling fraction versus inverse system 
size for Vi = 1.5 




FIG. 12. The superfluid density versus the interaction Vi 
(V2 = 0) and for different sizes. The systems are half filled 
and /3 = 6. 



we know that G is isotropic, and since we measured it 
along the lattice axis, we use these values to interpolate 
to lattice separations which we have not measured. For 
example, obtain G(-\/2), which corresponds to nnn, we 
interpolate the values of G{1) and G(2). The resulting 
Green function is Fourier transformed (two dimensional) 
to obtain the filling fraction. Figure 11a shows this fill- 
ing fraction for the fc^, = mode, psi^x — 0), for dif- 
ferent lattice sizes. We see a rapid transition around 
Vi = 2.1 from a phase where the condensate fraction is 
large (Vi < 2) to a phase where it is very small (Vi > 2). 
The inset of Fig. lib shows that the extrapolation of 
Psikx = 0) to the infinite system is linear in the in- 
verse system size. The extrapolated values are shown 
in Fig. lib exhibiting clearly the phase transition from 
non-zero to zero condensate density.c^ Figure 12 shows 
Ps versus Vi as directly measured from the winding num- 
ber. We see that the behaviour is exactly the same as 
for the condensate fraction, N{kx = 0), in Fig. 11. The 
transition of ps to zero goes hand in hand with that of 
the condensate fraction. These results for p^ and the 
zero mode filling fraction show that the transition from 
superfluid to (tt, 7r)-solid is continuous when V2 — 0. 

The question of what happens when this system is 
doped away from p — 0.5 is npl simple for large Vi and 
has been addressed elsewhere.E3 

So far we have discussed the two limits {Vi = 0, V2 7^ 0) 
and (Fi ^ 0,1^2 = 0). To construct the phase diagram 
in the {Vi,V2) plane we need to do simulations with 
both interactions non-zero. This was done in Ref. 12 
for Vi < 4.6, V2 < 3 where it was found that within these 
values of the interactions there was no direct (tt, 7r)-solid 
to (tt, 0)-solid transition although mean field calculatioas, 
predict this to happen starting at {Vi = 4, V2 = 2).t3 
Here we extend the simulations farther. Figure 13 shows 
Ps, S{Tr, tt), andS'(7r, 0) as a function of V2 for Vi = 6. We 
see that even at this large value of Vi there is no di- 
rect solid-solid transition. In addition, the sharpness of 
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V2 

FIG. 13. S'(7r,7r),S'(7r,0)andp, versus V2 for Vi = 6. There 
is no solid-solid transition even at this large value of Vi. 

the superfluid to (tt, 0)-solid transition suggests that it 
is first order while the relative smoothness of the super- 
fluid to (tt, 7r)-solid suggests that it is continuous. This 
is in agreement with Refs. 12 and 19. Preliminary re- 
sults indicate that even at Vi = 8 there still is no solid 
to solid transition. Apparently, the frustration caused 
by the competition between Vi and V2 prevents this and 
always results in a superfluid phase separating the two 
solid phases for any large but finite Vi^2- 

One of our goals in searching the phase diagram is to 
explore the possibility of a normal state at zero temper- 
ature. This is a particularly interesting question since 
such states have been found experimentally in high Tc 
superconducotrs. Explanations offered so far rely on the 
idea that the charge carriers in this normal state are the 
Cooper pairs, and thus Bosonic. However, the models 
which have been examined are for soft core Bosons at 
high densities. Our simulations for hard core Bosons have 
not shown any candidates for such a normal state. All 
phases we found so far are clearly solid, superfluid or 
solid as shown, for example, in Fig. 13. 
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VI. CONCLUSION 



1. General pcirtition function 



We have detailed the construction and testing of a 
Quantum Monte Carlo algorithm based on the exact du- 
ality transformation of the Bosonic Hubbard model. In 
particular we showed that this formulation yields a par- 
ticiilarly simple dual partition function in the very im- 
portant case of hardcore Bosons. We also showed how 
to measure various physical quantities. In particular, 
we showed how to measure the Green function, an im- 
portant quantity which is particularly difficult to obtain 
with other algorithms. To improve the statistics of this 
quantity we enhanced the efficiency of the algorithm by 
introducing a new reweighting procedure. 

We then used this algorithm to elaborate the structure 
of the various ground state phases of the Bosonic Hub- 
bard model. In particular, we showed and interpreted the 
behaviour of the Green function in the various superfluid, 
superolid and solid phases and calculated the filling frac- 
tion in the zero momentum mode. We found a filling 
fraction of about 33% at weak nn coupling, Vi = 1.5. 
This value is consistent with the pure hardcore value of 
Rcf. 21 (betwcxai 36% and 22%). However, the two re- 
sults were obtained at different Boson densities and at 
the moment it is not clear how our value depends on the 
density. This is in progress. 

We confirmed that the order of the phase transition 
from the superfluid to the solid phases depends on the 
long range order. At half filling, the SF-(7r, 7r)-solid tran- 
sition is second order while the SF-(7r, 0)-sohd transition 
is first order. 

The presence of a direct solid-solid transition from 
(tt, tt) to (tt, 0) is still an open question although the in- 
dications are that it is not present at any finite value of 
Vi^2- Another important issue is the possibility of a nor- 
mal conducting phase. Several different scenarios for the 
existence of normal conducting states in two dimensions 
have been offered but for softcore Bosons. So far, our ex- 
ploration of the phase diagram of this hardcore Bosonic 
model have given no hints of the presence of such a state. 
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APPENDIX A: DERIVATION OF THE DUAL 
PARTITION FUNCTION 

We present in some detail the integration of the orig- 
inal fields to obtain the dual expression for Z. Then we 

show, using hardcore bosons as an example, how to sim- 
plify this expression by explicitly solving constraints ap- 
pearing in it. 



To integrate out the bond variable 6 and the ampli- 
tude p in Eq. 29, we expand the exponentials of complex 
exponentials in Taylor series. This results in new integer 
variables /r,T, Pr)r and qi^l. I arises from the expansion 
of the second exponential in Eq. 29, p^'^ and g'^'^ from 
the i*'* hopping term. 



^=E n 



r,T r,T r,T 



,{~Pr,T-Pr,T) „-5V{nr,T) 



..(1) , Jl) 
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X Pr,T 

+Pfl2,r-1+Cfr}-1+---] 
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(2 

-rt^^^ -L ... - n _ A/^l 



_,-Bl r„(l)-„(l)+„(2) (2) _J2) 



X e 



where the sum runs over all the discrete indices we have 
introduced: I, p, q, J, N, and n. Integrating over 
6^, 0°and6° will give three constraints on the dual vari- 
ables enforced by Kronecker deltas. The constraint com- 
ing from the integration over 6^ ^, enforces the condition 



(A2) 



The total conserved current J^^ + N^^ is then just the 
current of bosons crossing the lattice. Negative time cur- 
rents are not allowed: Bosons cannot go backward in 
imaginary time. 

Using the two remaining constraints, one can see that 
the power of p^.r is 2(J°^_;^ -I- iV°). The integration 
over p then yields a (J^^ + N^)\ which compensates the 

l/n^^r! = l/{JrT + -^r)' already present in the expres- 
sion. 

The remaining p integrals are gaussian and can be per- 
formed easily. This leads to the dual formulation of Z: 



^=En 



(A3) 
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't-jT • Pr,T ■ HT,T • Pr,T ■ \lr,T ■ 



X r 1 + 



lr,r + lr,r+l + P^r-\,r + Pr^r+l + 



_ (2) 
1r-l,T 



where the sum runs over all the remaining integer in- 
deces I, p, q, and the conserved current J + N. 

It is impossible to construct an efficient Monte Carlo 
algorithm using this equation: There are many differ- 
ent fields, and many constraints. If one proposes a trial 
move, it will be rejected most of the time because it does 
not respect the conditions imposed by the S's or the fac- 
torials. 

Of course the role of the (5-functions is to allow only 
the acceptable configurations and give them the correct 
weight. It is possible to simplify greatly this partition 
function by solving these constraints, i.e. by finding the 
allowed configurations of currents and the corresponding 
weights. We now show how to do this for the case of 
hardcore bosons. 



2. Hardcore bosons 

To impose the hard core limit, we require the number 

of particles not be greater than 1 on any site. This is 
achieved simply by imposing that + = 0, 1 and 
I J^,- + = 0, 1 on all bonds. It is also understood 
that the currents must always be conserved. 

Now we need to find the Boltzmann weights of the al- 
lowed configurations. To do this we calculate the weights 
of representative Boson {i.e. current) configurations. 
The contributions of the interaction and chemical po- 
tential arc very simple to include. We have an cxp^Sfi) 
factor for each non-zero time current element, and an 
exp(— 5Vi) for each pair of near neighbour time currents. 

For example, consider configuration (a) in Fig. 13 
where the time current, J°^-|- Af°, is 1 along the directed 
path shown and zero everywhere else. This configura- 
tion satisfies all the constraints, and is uniquely specified 
by having all integer fields equal to zero except for the I 
field along the directed path where it is equal to 1. This 
is shown in the figure. So the Boltzmann weight of this 
configuration is simply 



[|e*^r(2) = (e''^)^- =e 



/3m 



(A4) 
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(a) (b) (c) 

FIG. 14. Different configurations of the conserved current, 
(a) without jumps, (b) with two single jumps, and (c) with 
one double jump and two single ones. The non-zero values of 
the integer dual variables are indicated. 

which agrees with intuitive expectations. 

Now consider path (b) in Fig. 13. Clearly, this is an 

allowed configuration and as before we determine, by ex- 
plicit calculation, the values of the integer fields that give 
it. These are shown in the figure. The Boltzmann weight 
of this configuration is therefore: 



(A5) 



The physical picture that emerges is simple: We always 
have the usual chemical potential factor and when the 
current line jumps from one site to another, we have a 
contribution of the hopping term. The I are non zero 
where there are time currents, the p where there is a 
jump to the right, and the q where there is a jump to the 
left. As a third example, to fix ideas, consider path (c) 
in Fig. 13. It has one double jump and two single ones. 
Therefore, we now have one non-zero p^'^^ , corresponding 
to the double jump to the right, and two q''^^ for the two 
single jumps to the left. The weight of this configuration 
is 



(A6) 



Thus, in the hardcore limit, the partition function sim- 
plifies drastically. Only one field survives: The conserved 
hard core current. The weight of a configuration can be 
easily calculated as follows 



Z=Y, V{N,J) 

{J,N} 



(A7) 
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where A'^^ is the number of jumps of length L, Nvi is the 
number of near-neighbor time current on the lattice, Nh 
is simply the number of bosons, and V is the Boltzmann 
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(a) (b) 

FIG. 15. Configuration with softcore bosons. (a) four 
bosons at tlie same site, (b) tliree bosons at tlie same site 
and one of them hopping to a neighboring site. 



weight. As mentioned in the main body of the paper, this 
partition function has systematic Trotter error of order 6. 
If one wants to construct an algorithm accurate to order 
(5^, one would have to keep track of all the second order 
terms all along the development of the kinetic term. 

One can generalize this approach to other kinds of sys- 
tems. For example if we look at the soft core configura- 
tions shown in Fig. 14(a), four bosons cross the lattice on 
the same site without hopping. The only non zero field 
is lr',T and we get the weight 



r(5)e 



„/3(4M-6yo) 



(A8) 



which is the expected result. Configuration (b) of Fig. 14 
with three bosons on the same site and one of them hop- 
ping to a neighboring site contributes 



3 e^''^ g-5yo(3iT 



-2.5^1 



(St)' 



(A9) 



The prefactor, 3, can be interpreted physically: Each 
boson can make this jump with equal probability and 
the sum must be taken. 
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